%% RunTravel is the base function of the TRAVEL model
% 
% input: cfg and gui, fromt he TRAVEL_MATLAB function or GUI
% ouput: saves plots and .mat files
%
% Defines all the parameters (in combination with the GUI's values), 
% sets the right triggers for the model
% performs the nonlinear optimizations,
% calls the functions of the submodels, 
% produces plots, 
% saves the data.
% 
function MainFunction(cfg,gui)
for n=[1 2]; if ishandle(n);close(n);end; end; % Closes some figures
set(gui.h.Play,'BackgroundColor',[0.2 0.9 0.2]);
pause(0.1); i=1; m=1;
data = cell(1,cfg.NoRuns);

%% Perform the required loops
[cfg,exvars,depvars,indepvars] = DefineParameters(cfg,gui);
cfg = MonteCarloSimulation(cfg,indepvars,1);

if isfield(cfg,'MonteCarloFields'); loopie = 1:max([cfg.LoopOverVariables... % Allows for looping over different variables independently in monte carlo
        *sum(cellfun(@(x) min([size(indepvars.(x),3),sum(indepvars.UseModes(1,:))]),cfg.MonteCarloFields)),1]);
else loopie=1; 
end

for variable = loopie
%% Loops over each variable
    cfg.variable = variable;
if cfg.LoopOverVariables
cfg.OutputFile = ...
    [cfg.MonteCarloFields{cfg.ParentFields(variable)} num2str(cfg.ModeFieldToUse(variable))];
end

if cfg.ParallelComputing==1 % Required for faster computations
if matlabpool('size') == 0; matlabpool(cfg.NumberofCPUCores); end
    parfor   Run = 1:cfg.NoRuns
        if ~cfg.NoOptimization % Avoids masses of data generation if no optimization is used
            [data{Run}.cfg data{Run}.indepvars data{Run}.depvars] = Iteration(cfg,gui,Run,exvars,depvars,indepvars);
        else
            [data{Run}.indepvars data{Run}.depvars] = Iteration(cfg,gui,Run,exvars,depvars,indepvars);
        end
    end

else
    for Run = 1:cfg.NoRuns
        if ~cfg.NoOptimization
            [data{Run}.cfg data{Run}.indepvars data{Run}.depvars] = Iteration(cfg,gui,Run,exvars,depvars,indepvars);
        else
            [data{Run}.indepvars data{Run}.depvars] = Iteration(cfg,gui,Run,exvars,depvars,indepvars);
        end
    end
end

%% Save the data to file
[junk1 junk2] = mkdir(['Output/' cfg.OutputFolder]);
if ~cfg.NoOptimization
    save( ['Output/' cfg.OutputFolder '/' datestr(now,'yyyy-mm-dd_HH-MM-SS' ) '-'  cfg.OutputFile '.mat'] ,'data');
else
    data{1}.cfg = cfg;
    save( ['Output/' cfg.OutputFolder '/' datestr(now,'yyyy-mm-dd_HH-MM-SS' ) '-'  cfg.OutputFile '.mat'] ,'data');
end
set(gui.h.Play,'BackgroundColor',[0.8 0.8 0.8])
disp(['Saved to Output/' cfg.OutputFolder '/' datestr(now,'yyyy-mm-dd_HH-MM-SS' ) '-'  cfg.OutputFile '.mat']);
end

disp('Done');              

end

function varargout = Iteration(cfg,gui,Run,varargin)
%% Set configuration parameters and preloads variables

if nargin > 3
    exvars    = varargin{1};
    depvars   = varargin{2};
    indepvars = varargin{3};
end

if ~exist('depvars','var')
    [cfg,exvars,depvars,indepvars] = DefineParameters(cfg,gui);
end

cfg.Run = Run; 

[cfg,exvars,depvars,indepvars,cfg.Regions] = TransformModesandRegions(cfg,exvars,depvars,indepvars,cfg.RegionSet);
[cfg,exvars,depvars,indepvars]             = TransformModesandRegions(cfg,exvars,depvars,indepvars,cfg.ModeSet);
if cfg.AllRegions == 1; cfg.CalculatedRegions = cfg.Regions; end

    for r = cfg.CalculatedRegions
        %% Randomly vary inputs (optional)
        if cfg.ParallelComputing==0; set(cfg.h.Info,'String',{'Region / Time / Run: ',[num2str(r) ' / ' num2str(1) ' / ' num2str(cfg.Run)]}); end
        [cfg,indepvars] = MonteCarloSimulation(cfg,indepvars,r,cfg.variable*cfg.LoopOverVariables);

        %% Calibrate
        if cfg.NoOptimization % In case of no optimization
        elseif cfg.MonteCarloSimulation && cfg.VanRuijvenMonteCarlo
            disp (['Running ' cfg.Model ', region: ' num2str(r) ', run: ' num2str(cfg.Run) ' Optimization over Monte Carlo set']);
            [cfg,indepvars, depvars] = OptimizeAfterMonteCarlo(cfg,exvars,depvars,indepvars,r);
        else
            disp (['Running ' cfg.Model ', region: ' num2str(r) ', run: ' num2str(cfg.Run) ]);
            [cfg,indepvars, depvars] = Calibration(cfg,exvars,depvars,indepvars,r);
        end

        %% Project
        depvars = Projection(cfg,exvars,depvars,indepvars,r);
        %% Plot
        GeneralPlotting(cfg.plotmode,'save',cfg,exvars,depvars,indepvars,r,[]);
    end
    
if nargout ==3 || ~isfield(cfg,'MonteCarloFields')
    varargout = {cfg, indepvars, depvars};
else
    umcfields = unique(cfg.MonteCarloFields);
    for nn = 1:length(umcfields)
        indepout.(umcfields{nn}) = indepvars.(umcfields{nn});
    end
    
    relfields = {'TravelDemand','ModeSplit','errors'};
    for nn = 1:length(relfields)
        depout.(relfields{nn})   = depvars.(relfields{nn});
    end
    varargout  = {indepout,depout};
end
end


function [cfg, indepvars] = MonteCarloSimulation(cfg,indepvars,r,varargin)
%% Runs the calibration procedure as deinfed by Bas van Ruijven
% Ranges can be defined in the cfg.inputcell part below
modes    = indepvars.UseModes(r,:,:);
nummodes = find(modes);
if (nargin < 4 || varargin{1}==0)&& isfield(cfg,'SelectionResults') 
    varargin{1} = 1:length(cfg.ModeFieldToUse);
end


if cfg.MonteCarloSimulation ==1
if ~isfield(cfg,'SelectionResults')
    %% Run it the first time to set some stuff up
        switch cfg.Model
            case 'ModelGcam'
                cfg.inputcell = {'CostDistrParameter', 'PriceElasticity','IncomeElasticity', 'Sigma','ShareWeight'};
                ranges    = {[0.9 0.2], [0.55 0.15], [1 0.2], [0 0],[1 0.2]};
                
            case {'ModelTravelOptimizedTW','ModelTravelLinearTW'}
                cfg.inputcell = {'Lambda','TTB','TMBAdded','ComfElasticity','Inertia','PreferenceFactorSlope','PreferenceFactorStart','TimeWeight'};
                ranges    = {[3.2 0.5], [400 50], [0 0.01],[0.5 0.1], [1 0.66], [1 0.05],[1 0.8],[0 0.1]};
            case {'ModelPoles'}
                cfg.inputcell = {'IncomeElasticityER','TimeTrend','ERSaturation','PriceElasticity','IncomeElasticity'};
                ranges        = {[1 0.2],[0.005 0.0001],[1 0.1],[1 0.1],[1 0.1]};
                
            case {'ModelSimple'}
                cfg.inputcell = {'PriceElasticity','IncomeElasticity'};
                ranges        = {[1 0.1],[1 0.1]};
                
            case 'test' 
               cfg.inputcell = 0;
        end
    
   %% Generate requires GUI elements 
   
        cfg.SelectionResults = SelectionMenu(cfg.inputcell);
        regions              = str2double(cfg.RegionSet(end-1:end))+1;
        cfg.MonteCarloFields = cfg.inputcell(cfg.SelectionResults);
        cfg.ranges           = ranges(cfg.SelectionResults);
    
        configs.mcnumel   = cellfun(@(x) size(indepvars.(x),3),cfg.MonteCarloFields);
        configs.mcnumel(configs.mcnumel==length(cfg.modes)) = sum(modes); 
        configs.variables = {}; configs.variablesnum = {}; configs.tranges   = [];
        cfg.ModeFieldToUse= [];
        cfg.ParentFields  = [];
        
        z=1;m=1;mi=1;i=1;
        for n = 1:sum(configs.mcnumel)
            if configs.mcnumel(min([i,mi]))>1 && m<configs.mcnumel(min([i,mi]))
                z = z+0.99;
            elseif configs.mcnumel(min([i,mi]))>1
                z=1.99; mi=mi+1;
            else
                z=1.99; mi=mi+0.99;
            end
            i=floor(mi);
            m=floor(z);
            configs.variables(end+1)    = cfg.MonteCarloFields(i);
            configs.tranges(end+1,1:2)  = cfg.ranges{i};
            if configs.mcnumel(i) == 1 && size(indepvars.(cfg.MonteCarloFields{i}),3)==1;
                cfg.ModeFieldToUse(end+1)   = 1;
                txt=''; 
            else
                txt=num2str(nummodes(m)); cfg.ModeFieldToUse(end+1)   = nummodes(m);
            end
            configs.variablesnum{n}     = [cfg.MonteCarloFields{i},  txt];
            
            cfg.ParentFields(end+1)     = i;
        end
               
        configs.ranges    = NumberEntryMenu(configs.variablesnum , num2cell(configs.tranges,2));
        
        cfg.validranges   = repmat(configs.ranges,[1,max(cfg.CalculatedRegions)]);

end
    
if ~isfield(cfg,'MonteCarloDesign')
    %% Create the design for the monte carlo experiment
    if sum(cfg.SelectionResults)>1
        cfg.MonteCarloDesign = ccdesign(size(cfg.validranges,1),'type','inscribed','fraction',0,'center',2);
    else
        cfg.MonteCarloDesign =[];
    end
    
    if size(cfg.MonteCarloDesign,1)>cfg.NoRuns; 
        cfg.MonteCarloDesign = 2*(rand(cfg.NoRuns,(size(cfg.validranges,1)))-0.5);
        disp(['DCC designed with ' num2str(size(cfg.MonteCarloDesign,1)) ' random variables' ]);   
    else 
        disp(['DCC designed with ' num2str(size(cfg.MonteCarloDesign,1)) ' CCC and ' num2str(cfg.NoRuns - size(cfg.MonteCarloDesign,1)) ' random variables' ]);   
        cfg.MonteCarloDesign(end+1:cfg.NoRuns,:) = 2*(rand(cfg.NoRuns - size(cfg.MonteCarloDesign,1),(size(cfg.validranges,1)))-0.5);
    end
    
else
    %% Transform the input variables
    temp.info           = {};
            
    for m = varargin{1};
        indepvars.(cfg.MonteCarloFields{cfg.ParentFields(m)})(r,:,cfg.ModeFieldToUse(m)) =...
            cfg.validranges{m,r}(1) + cfg.validranges{m,r}(2)*cfg.MonteCarloDesign(cfg.Run,m);
        eval(['temp.info(end+1) = {[' '''; '''  ' cfg.MonteCarloFields{cfg.ParentFields(m)} num2str(cfg.ModeFieldToUse(m)) '...
            ''' = ''' ' num2str(mean(indepvars.' cfg.MonteCarloFields{cfg.ParentFields(m)} '(r,:,cfg.ModeFieldToUse(m))),2)]};'] );
    end
    cfg.plotinfo{r} = char(temp.info);
    cfg.plotinfot   = cfg.plotinfo{r}';
    
    if ~cfg.NoOptimization
         disp (['Monte Carlo ' cfg.plotinfot(:)']); 
    else
         if mod(cfg.Run,100) == 0; disp(['Run number: ' num2str(cfg.Run)]); end
    end
end
end
end